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Abstract —  Continuity  equation  for  wave  modeling  is  still  being  developed.  There  are  quite  a  lot  of  versions  of  this 
equation.  This  research  formulates  continuity  equation  in  a  simple  form  to  simplify  its  numerical  and  analytical 
solution. 

The  formulation  of  the  continuity  equation  is  done  by  performing  mass  conservation  law  in  a  water  column  with  free 
surface  and  by  performing  weighted  total  acceleration.  Then,  the  continuity  equation  is  performed  along  with  the 
surface  momentum  equation  and  completed  numerically  to  modeling  one-dimensional  wave  dynamism.  The  equation 
is  capable  of  modeling  shoalingand  breaking. 
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I.  INTRODUCTION 

Time  series  water  wave  equation  is  generally  called 
Boussinesq  type  equation.There  are  quite  a  lot  of  versions  of 
Boussinesq  equation,  either  its  continuity  equation  or  water 
surface  equation  as  well  as  its  momentum  equation.  Those 
equations  generally  consist  of  the  second  or  higher 
differential  elements  which  quite  complicate  the  solution 
both  analytically  and  numerically.  Some  researcher  who  have 
developed  Boussinesq  equation  are  among  othersBoussinesq, 
J.  (1871),  Dingeman,  M.W.  (1997),  Ham;L„  Madsen,  P.A., 
Peregrin,  D.H.  (1993),  Johnson,  R.S.  (1997),  Kirby,  J.T. 
(2003),  Peregrine,  D.H.  (1967),  Peregrine,  D.H.  (1972)  and 
many  more. 

Governing  equations  in  this  research  are  water  surface 
equation  and  surface  momentum  equation,  both  of  which  use 
particle  velocity  at  the  surface  as  the  variable  and  both  are  in 
the  form  of  time  and  space  differential  equation  in  simple 
form.  Water  surface  equation  is  formulated  based  on  mass 
conservation  law  and  by  performing  weighted  total 
acceleration  at  the  kinematic  free  surface  boundary 
condition.  The  momentum  equation  is  obtained  by 
performing  weighted  total  acceleration  at  the  Euler 
momentum  equation.  The  integration  with  water  depth  from 
this  momentum  equation  produces  surface  momentum 
equation  with  particle  velocity  variable  at  the  surface. 


Both  governing  equations  are  done  using  numerical  method 
where  spatial  differential  is  done  using  finite  difference 
method,  whereas  time  differential  is  done  using  corrector 
predictor  method. 


H.  TOTAL  DERIVATIVE  EQUATION 

Hutahaean  (2019a)  developed  weighted  total  acceleration 
equation  at  water  particle  in  horizontal  direction,  i.e.. 


Du  3u  ,  3u  3u  ... 

—  =  V - h  U - hW —  . (1) 

uis  particle  velocity  in  horizontal-xdirection  andwis  particle 
velocity  in  vertical-z  direction.  Weighted  total  acceleration, 
was  actually  formulated  for  the  function  /  =  f(x,t). 
However,  in  this  research  it  is  performed  at  f  —  f(x,z,t), 
because  the  wave  being  discussed  is  a  wave  moving  to 
horizontal-xdirection  and  vertical  z  dimension  is  eliminated 
with  the  integration  process,  so  the  equation  becomes  a 
function  of  /  =  fix,  t). 


The  changes  in  total  water  surface  elevation  is 

Bl  =  ya-1  +  U . (2) 

dt  'at  'i  ax  w 

7]  —  rj{x,  t)is  water  surface  elevation  against  still  water  level 
(Fig.  1).  In  (1)  and  (2)  there  is  time  coefficient  or  time  scale 
at  time  differential  i.e.  ywith  a  value  of  2.87-3.14  in 
Hutahaean  (2019a,b).  The  value  of  y  is  very  much 
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determined  by  basic  equation  in  which  the  total  acceleration 
was  performed.  In  this  research  the  correspond  y  value  of 
2.00  is  found. 


Based  on  (2),  then  kinematic  free  surface  boundary  condition 
that  was  formulated  from  total  derivative  equation  of  the 
changes  of  a  surface  (Dean  (1991)),  becomes 


w„ 


a>7 


ai?  . 

y — h  u„ 
'at  ''ax 


Vv  ~  r  3,  1  “r;  . (3) 

uv is  the  velocity  of  horizontal -^direction  at  the  surface. 


III.  CONTINUITY  EQUATION 

3.1.  The  formulation  of  continuity  equation 
Continuity  equation  or  water  wave  surface  equation  will  be 
formulated  in  a  water  column  (Fig.l.)  with  free  water 
surface,  where  as  a  result  of  an  input-output  in  a  water 
column  in  a  very  small  time  interval  ySt,  a  change  in  water 
surface  elevation  of  Sr]  occurs  so  there  is  also  a  change  in  the 
per  width  unit  volume  of  Sr/Sx.  For  a  very  small 
GxwhereGx  =  dx 
Sm  —  j Sr/dx  . (4) 

The  change  in  water  mass  from  input-output  process  (Fig.  1.) 
is, 

Sm  =  p(u  —  (u  +  Su^SzSt 

+  (w  —  (w  +  Sw))SxySt 
rSu  5w\ 

Sm  =  —  p  —  +  ——  SxSzySt 
\Sx  SzJ 


Total  change  in  water  mass  in  the  water  column  at  very  small 
SxandSz, 

Sm  =  -p  dz  dxySt  . (5) 

his  water  depth  against  still  water  level,  T]  —  r](x,t) is  the 
water  surface  elevation  also  against  water  level. .For 
incompressible  fluid,  Smin  (4)  is  the  same  asb'min  (5), 


p  f77  /du  3w\ 

-Sijdx  =  -p  j^—  +  —]dzdxySt 

Both  sides  are  divided  by  p,  dxandySt, 


For  a  very  small  fit, 

-LU_f 

2 ydt  ]_h  V3x  3z/ 

Substitute  (2)  to  the  left  side  of  the  equation 

1  a p  uv  ap  r 77  /a u  aw\ 

2  at  +  2y  3x  J_h  Va^  +  a  z)^Z 

The  integration  of  the  second  term  of  the  right  side  is  done 
and  substituted  kinematic  free  surface  bondary  condition  and 
kinematic  bottom  boundary  condition. 


The  integration  of  the  first  term  right  side  is  performed  with 
Leibniz  integration  (Protter  (1985)), 


f 

J  a 


a  /  a 

—  dz  =  — 
3x  3x 


Obtain, 

^  3u  a 

—  dz  -  — 
dx  dx 


fP  3 B  da 

f 

J-h 


37? 


XL  d.Z  XLyt  „  XL_u 


L 

(6),becomes 

+  udZ-^ 

V  2/ at  axJ-h  2v  ax 


3  h 
dx 


■(7) 


Integration  in  the  right  side  in  (7)  is  done  using  velocity 
equation  from  Dean  (1991)  and  the  result  of  integration  is 
expressed  as  a  function  of  surface  horizontal  velocity  iirj  in 
order  to  correspond  with  momentum  equation  that  produces 
surface  velocity u^.  Velocity  potential  equation  as  the  result 
of  Laplace  equation  solution  (Dean  (1991))  is 
0(x,  z,  t)  =  GcoskxcoshkQi  +  z)sinot  ....(8) 

Gwave  constant,  k  wave  number  ander  angular  frequency. 
Particle  velocity  in  horizontal-x  direction  is 
3  cP 

u(x,z,t )  =  —  — - 
dx 

=  Gksinkxcoshk(h  +  z)sin<jt . (9) 


T  T  .  coshk(h+z ) 

Using  (9),  U  =  - ttt: - -Uyy 


coshk(h+ri')  ^ 


then  integration  in  (7) 


becomes. 


a  n  a  f 

^Ludz  =  5*1 


77  coshk(h  +  z) 


_h  ux j_hcoshk{h  +  p)  71 
Completing  the  integration  will  obtain 


i i-H  dz 
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3  n  a  / u„tanhk(h  +  rj) 

—  I  udz  —  —  I - 


3 xj_h  3x  y 

From  the  wave-number  conservation  equation  (Hutahaean 
(2019a)),  tanhk(h  +  g)  —  tanhk0(h0  +  g0)  =  1,  where/c0is 
wave  number  in  deep  water,  h0 is  deep  water  depth  andt70is 
water  surface  elevation  in  deep  water,  can  have  a  value  of 
yor  others,  A0is  wave  amplitude  in  deep  water.  Therefore, 
the  result  of  the  integration  becomes, 

a  r11  a  /ii„N 

sj_,“dz  =  5t(  T) 

Substitute  the  result  of  integration  to  (7), 


V  2/  at  a x\k  J 


U77  a?7 
2  y  a% 


.(10) 


Equation  (10)  is  a  continuity  equation  that  will  be  used  in  this 
researchor  water  wave  surface  equation  in  the  form  of 
differential  equation.  In  (10),  there  is  wave  number 
^parameter  that  should  be  known,  and  some  other 
characteristics  that  should  also  be  known,  among  other  is 
deep  water  depthd0,  i.e.  maximumwater  depth  if  the  equation 
was  done  in  water  depth  c/which  is  bigger  than  d0,  so  the 
calculation  is  done  using  cl0. Next  is  wave  amplitude 
rnaximum/lmax,  i.e.  maximum  amplitude  in  a  wave  period 
that  can  be  inputted  to  the  model. 


3.2.  The  calculation  oT4maxandd0. 

It’s  been  known  that  there  is  a  relation  between  water  depth 
c/and  wave  number  k,  then  the  calculation  will  be  easier  if  in 
(10)  wave  number  kis  substituted  with  water  depth  d. 
Whereas  the  equation  for  wave  number  in  deep  water  k0can 
be  calculated  using  the  following  equation,  the  formulation 
of  an  equation  outside  the  scope  of  this  research,  will  be 
written  in  the  next  paper. 

y(y +  l)°2  =  gk0(l- k0A0)  . (11) 

A0is  wave  amplitude  which  is  an  input,  k0 is  deep  water  wave 

2n 

number,  o  is  angular  frequency,  a  —  — ,  Tis  wave  period. 
k0in(ll)  can  be  calculated  using  simple  calculation,  i.e. 
finding  the  root  of  the  quadratic  equation.  (11)  can  be 

completed  if  determinant  D  =  g2  —  AgAy  (y  +  ^  <r2is 
bigger  than  or  the  same  as  zero. In  case  of  D  —  0, obtains 

^0 ,max  ~  — 7  FT~r  . (12) 

4y(y+jJo- 

In  deep  water  tanhk0  ( d0  +  y)  =  1  applies.  Assuming  that 
wave  amplitude  is  much  smaller  than  deep  water  depth,  or 


—  «  1,  A0  deep  water  wave  amplitude  andd0  deep  water 

2d0 

depth,  then  the  following  relation  applies 

/  A>\  (  A0  \ 

tanhk0  \  d0  +  —J  —  tanhk0d0  (^1  +  J 

= tanhk0d0  =  1 

k0is  wave  number  in  deep  water  depth  d0.  As  deep  water  the 
following  criteria  is  used 
k0d0  =  1.77 r  . (13) 

where tanhl.Vn  —  0.999954  «  1,  the  uses  of  this 
1.7  7rvalue  is  also  based  on  the  review  of  the  produced 
breaker  depth.  k0 was  obtained  from  (11),  therefored0can  be 
calculated  using  (13). 


Bases  on  wave  number  conservation  equation  (Hutahaean 
(2019a)),  the  relation  between  wave  number  kd  in  a  depth  d, 
with  wave  number  in  deep  water  is, 

kd(h  +  r/)  —  k0(d0  +  r]0).  Assuming  that  —  «  1,  then 

do 


kd(h  +  rj)  =  k0d0  = 


kd 


l.ln 

(h+J?) 


(14) 


1.7 n  .Or, 


3.3.  The  Final  Water  Wave  Surface  Equation 
By  substituting  (14)  to  (10),  the  final  water  wave  equation 
was  obtained  with  water  depth  das  its  parameter,  i.e 
0.  ,¥1-  a  (uyV+m  uv  ar, 

V  2/  a t  3x  V  1.7JT  J  2y3x  . ^  ' 

Therefore.  There  is  no  need  to  calculate  wave  number  k. 

An  example  of  the  calculation  of  deep  water  wavelength 

2n 

L0  —  —  ,  deep  water  depth  d0andA0max where  y  —  2.000is 
used  is  presented  in  Table  (1)  below. 


Table. 1:  The  value  of  A0  max,  d0andL0 


T 

A 

n0,max 

d0 

Lo 

d0 

(sec.) 

(m) 

(m) 

(m) 

^0 

6 

0,447 

4,929 

5,798 

0,85 

7 

0,608 

6,708 

7,892 

0,85 

8 

0,794 

8,762 

10,308 

0,85 

9 

1,005 

11,09 

13,047 

0,85 

10 

1,241 

13,691 

16,107 

0,85 

11 

1,502 

16,566 

19,489 

0,85 

12 

1,787 

19,715 

23,194 

0,85 

13 

2,098 

23,137 

27,221 

0,85 

14 

2,433 

26,834 

31,569 

0,85 

15 

2,793 

30,804 

36,24 

0,85 
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..(17) 

••(18) 


IV.  MOMENTUM  EQUATION 

Weighted  total  acceleration  equation  is  done  in  Euler 
momentum  equation  in  horizontal-xdirection  and  vertical- 
zdirection  consecutively  (Anderson  (1995)), 

3u  3u  3u  1  3p 

V - b  U - hw —  = - 

at  ax  az  p  ax 

aw  ,  aw  ,  aw  1  a p 

'  at  ax  az  p  az  a 

(18)  is  written  as  an  equation  for  pandthe  nature  of  irrotional 

flow  was  performed  i.e.^f  =  |p,  the  equation  is  integrated  to 
vertical-z  axis,  surface  dynamic  boundarycondition  is 
performed,  i.e.  p);  =  0,  and  differentiated  against  horizontal- 
x  axis. 

lap  a  r?'3w  .  i  a 

pi 

(19) 


i  aP  a  paw  1  a  2, 

-y—\  — —  dz  + -— (u2  +  w2) 
p  3x  1  3x  J.  3t  2  3*v  v 


-(u2  +  w2)  +  g- 
2  ax  K  J  a  ax 

In  (17)  the  nature  of  irrotional  flow  was  performed,  i.e.  ^  = 
and  substitute  (19)  to  the  right  side  of  the  equation. 


au  a  pi  aw 

rai  =  ~rail  -aidz 

1  a  (un  +  WP  -  u  'tz . (20) 


2  ax 


The  solution  of  ^  dzis  done  using  velocity  potential 

(8),  where  the  particle  velocity  in  horizontal  direction  is  in 
equation  (9).  Particle  velocity  in  vertical-z,  is 

a$ 

=  —Gksinhk(h  +  z)coskxsinot 


3  z 


=  —Gksinhk(h  +  z)<jcoskxcos<jt 


....(21) 
aw 
at 


.(22) 


(22)  is  integrated  against  time  t 


l 


11  a  w 

— — dz  = 

at 


-dz  — 


—G{coshk(h  +  p)  —  coshk(h  +  z)) 
ocoskxcosot 

Differentiated  against  horizontal-*  axis 

a  paw 
a x  Jz  at 

Gk(coshk(h  +  p)  —  coshkQi  +  z)) 
osinkxcosot 

Equation  (9)  is  differentiated  against  time  t,  —  = 

at 

Gksinkxcoshk(h  +  z)acosat ,  which  shows  that  this  form  is 

in  —  f’1  —  dz,  so  the  following  relation  is  obtained 
ax  Jz  at  b 


a  r11 3w  3u„  a u 

—  - dz  =  — - - 

a x  Jz  at  at  at 

Substitute  this  equation  to  (20), 


—  =  -  (u2  +  w2)  +  g— )-  . (23) 

at  V2  ax v  '  a  a x)  y  v  ' 

(23)  is  surface  momentum  equation  that  produces  surface 


velocityu);.  By  completing  “jy dzwith  the  method 

above,  then  in  the  momentum  equation  there  is  an  influence 
of  continuity  equation  or  momentum  equation  that  was 
produced  and  controlled  by  continuity  equation.  Another 
control  by  continuity  equation  on  momentum  equation  is  on 


variablew„  ,  i.e.w„  =  v  —  +  —where— is  obtained  from 

v  v  '  at  71  ax  at 

continuity  equation.  Therefore,  momentum  equation  (20)  is 

controlled  by  water  surface  equation. 


V.  RESULT  OF  THE  MODEL 

5.1.  Numerical  Solution 

In  this  research,  water  surface  equation  and  momentum 
equation  are  done  with  finite  difference  method  for  spatial 
differential,  whereas  time  differential  is  done  using  predictor- 
corrector  method  based  on  Newton-Cote  numerical 
integration  (Abramowitz  (1972)).  Whereas  the  predictor- 
corrector  method  is  as  follows.  As  an  example  water  surface 
equation  (15)  will  be  done.  The  water  surface  equation  can 
be  written  in  the  form  of. 


The  equation  is  integrated  against  time  from  t  —  t  —  btuntil 
t  —  t  +  St,  where  the  integration  of  the  right  side  of  the 
equation  is  done  with  Newton-Cote  numerical  integration 
with  3  (three)  integration  points, 

rt+St  rt+St 

I  3  p  =  I  F(t)dt 
Jt-St  Jt—St 

rjt+St  =  rf-st  +  St  pt-st  +  j_Ft  +  ipt+stj 

....(25) 

Ft+St is  unknown  number,  therefore  it  needs  to  be  predicted 
using  Taylor  series  and  finite  difference  method,  where  the 
step  is  called  predictor  step,  i.e. 

Ft+St  =  Fl  +  St  ....(26) 

Substitute  (26)  to  (25),  the  value  of  pt+5tprediction  can  be 
calculated.  With  similar  way,  momentum  equation  is  done, 
and  prediction  is  obtained.  With  those  prediction 

values, Ft+<Stcan  be  calculated  with  (24),  and  (25)  is  done. 
This  step  is  called  corrector  step.  This  corrector  step  is  also 
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done  on  momentum  equation  interchangeably  with  the  water 
surface  equation,  and  repeated  until  a  convergence  is  reached 
where  |  -  Votd1  |  <  eand|u^n5/w  -  \  <  swheresis 

a  very  small  positive  number  as  an  iteration  accuracy  criteria. 
Generally,  a  convergence  is  reached  with  5  iterations. 

5.2.  The  Result  of  Model  Execution 
a.  In  constant  depth  of  1 1 .0  m 

In  this  section  the  model  is  done  in  a  channel  with  a  constant 
water  depth  of  d  —  11.00  m,  with  wave  period  of  8  sec., 
wave  amplitude  A0  —  0.794  m,  where  actually  that  does  not 
mean  that  wave  height  is  twice  that  of  wave  amplitude, 
Hutahaean(2019  c).Deep  water  depth  for  this  wave  is  d0  — 
8.762m.  In  the  case  that  dis  bigger  thand0then  the 
calculation  of(d  +  ?7)in  (15), (d0  +  rf)  is  used.  The  model  is 
done  using  two  boundary  conditions,  i.e.  closed-end 
boundary  condition  where  horizontal  velocityu  =  0,  whereas 
in  the  opened-endthe  model  was  given  an  input,  i.e. 
sinusoidal  wave?70  =  A0sirujt.  The  input  is  done  only  for  1 
time  wave  period. 

The  result  of  the  execution  for  8,  24,  and  40  sec.  is  presented 
in  (Fig. 2.).  In  the  execution  for  8  sec.,  the  wave  profile  is  still 
in  the  form  of  sinusoidal,  but  the  wave  trough  elevation  is 
smaller  than  the  elevation.  In  the  execution  for  24  sec,  the 
formed  wave  trough  is  getting  smaller  and  farther  away, 
similarly  with  execution  for  40  sec,  the  wave  trough  is 
getting  smaller  and  farther  away  where  water  ripple  is 
formed  and  the  form  of  the  main  wave  is  a  perfect  cnoidal 
wave  or  more  accurately  it  is  called  solitary  wave. 

As  a  conclusion  of  the  model  execution  in  this  constant  water 
depth  is  that  in  the  deep  water,  the  equation  used  produced 
perfect  cnoidaltype  wave  or  also  can  be  called  as  solitary 
wavetype,  even  though  the  input  of  sinusoidal  wave,  the 
wave  trough  part  disappears. 


t=8.0  sec.  t=24.0sec.  t=40.0  sec. 

Fig.2.Wave  profile  after  the  execution  for  8,24  and  40  sec. 


b.  In  a  changing  depth 

With  the  phenomenon  of  the  evolution  of  sinusoidal  wave 
into  cnoidal  wave,  in  the  model  execution  at  the  an  uneven 
bottom,  before  the  wave  enters  the  water  with  slopping 
bottom,  the  wave  is  given  evolution  zone,  i.e.  in  front  of  the 
water  in  the  form  of  water  with  constant  depth. 


Evolution  zone 


Fig. 3.  Sea  bed  for  model  execution  at  uneven  bottom. 


The  calculation  is  done  with  evolution  zone  length  of  100  m 
(Fig. 3)  with  constant  water  depth  old  =  11.0  m,  then  the 
water  depth  changes  until  the  depth  of  1 .0  m,  with  a  distance 
of  200  m,  with  tangent  of  bottom  slope  i.e.  —  =  —  =  0.05. 
The  wave  used  here  is  wave  with  wave  period  8  sec.,  wave 
amplitude  0.794  m,  with  the  result  of  calculation  shown  in 
Fig.4.  and  Fig. 5. Coming  out  of  the  evolution  zone,  shoaling 
occurs  followed  by  breaking,  with  a  breaker  height  Hb  — 
1.546  m,  at  breaker  depth  hb  —  1.969  m,  where  —  =  0.785. 

hb 

This  condition  is  obtained  by  multiplying  the  second  term  of 
the  water  wave  surface  equation  (15)  with  a  factor  of  2.5,  so 
that  (15)  becomes, 

(  |  11  a>7  _  a  fVd+qA  2-5uv  ay 
v  2/  at  ax  v  i.7jt  J  2 y  ax' 

This  coefficient  2.5  is  obtained  by  experimentation  in  order 
to  obtain  —approximates  0.80.  Coefficient  2.52  can  also  be 

hb 

used  where  ^  =  0.81was  obtained  but  the  equation  becomes 

unstable  after  the  breaking.  Therefore,  further  research  is  still 
needed  on  water  wave  surface  equation  as  well  as  momentum 
equation  that  was  used. 
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Fig  4.  Wave  profile  (p)  and  wave  crest  elevation(rjmax ) 
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Fig  5.  Wave  profile  (tj)  after  breaking 

VI.  CONCLUSION 

As  has  been  shown  that  model  can  produce  two  main 
phenomena  that  occur  at  the  water  wave  on  its  way  to 
shallow  water,  i.e.  shoaling  and  breaking.  At  the  deep  water, 
at  a  constant  depth  the  profile  of  perfect  cnoidal  wave  is 
formed  which  is  also  called  solitary  wave.  Behind  the  main 
wave,  wave  ripple  is  formed  which  is  also  known  as  undular 
wave.  Therefore,  it  can  be  said  that  the  equation  that  was 
produced  in  this  research  can  model  several  phenomena  at 
water  wave  found  in  the  nature. 

Further  research  that  needs  to  be  done  is  to  study  the 
phenomenon  at  the  equation  by  producing  analytical 
solution.  Considering  the  simple  form  of  the  equation,  the 
analytical  solution  of  the  water  wave  surface  equation  can  be 
obtained  easily,  i.e.  using  velocity  potential  equation  from 
Laplace  solution  equation.  By  studying  analytical  solution,  it 
is  expected  that  an  explanation  will  be  obtained  on  the 
appearance  of  coefficient  2.5  at  the  second  term  of  the  water 
wave  surface  equation. 
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